#include "IsoBandUtil.h"

double earthRadius = 63710088e-1;
double FACTOR = (earthRadius * earthRadius) / 2;
double PI_OVER_180 = M_PI / 180;
void ringArea(const std::vector<IsoBandPoint> &coords, double &outArea, std::vector<double> &outBbox)
{
    outArea = 0.0f;
    if (coords.size() <= 3)
    {
        return;
    }
    outBbox.resize(4);

    double west = 180.0;
    double south = 90.0;
    double east = -180.0;
    double north = -180.0;

    size_t coordsLength = coords.size() - 1;
    double total = 0.0;
    size_t i = 0;
    while (i < coordsLength)
    {
        const IsoBandPoint &lower = coords[i];
        const IsoBandPoint &middle = coords[i + 1 == coordsLength ? 0 : i + 1];
        const IsoBandPoint &upper = coords[i + 2 >= coordsLength ? (i + 2) % coordsLength : i + 2];
        double lowerX = lower.x * PI_OVER_180;
        double middleY = middle.y * PI_OVER_180;
        double upperX = upper.x * PI_OVER_180;
        total += (upperX - lowerX) * sin(middleY);
        i++;

        west = std::min(lower.x, west);
        south = std::min(lower.y, south);
        east = std::max(lower.x, east);
        north = std::max(lower.y, north);
    }
    outBbox[0] = west;
    outBbox[1] = south;
    outBbox[2] = east;
    outBbox[3] = north;
    outArea = abs(total * FACTOR);
}

double linear(double a, double b, double v)
{
    if (a < b)
        return (v - a) / (b - a);
    return (a - v) / (a - b);
}

double interpolate(double a, double b, double v0, double v1)
{
    double tmp;

    if (v0 > v1)
    {
        tmp = v0;
        v0 = v1;
        v1 = tmp;
    }

    if (a < b)
    {
        if (a < v0)
            return (v0 - a) / (b - a);
        else
            return (v1 - a) / (b - a);
    }
    else if (a > v1)
    {
        return (a - v1) / (a - b);
    }

    return (a - v0) / (a - b);
}

double interpolate_a(double a, double b, double minV, double maxV)
{
    if (a < b)
        return (minV - a) / (b - a);
    return (a - maxV) / (a - b);
}

double interpolate_b(double a, double b, double minV, double maxV)
{
    if (a < b)
        return (maxV - a) / (b - a);
    return (a - minV) / (a - b);
}


// 创建指定长度的向量
Vec4 vec4()
{
    return Vec4(4);
}

Vec8 vec8()
{
    return Vec8(8);
}

Vec12 vec12()
{
    return Vec12(12);
}

Vec16 vec16()
{
    return Vec16(16);
}

// 计算向量元素和
double estimate(int elen, const Vec4 &e)
{
    double Q = e[0];
    for (int i = 1; i < elen; i++)
    {
        Q += e[i];
    }
    return Q;
}

// 求和函数，处理精度问题
int sum(int elen, const Vec4 &e, int flen, const Vec4 &f, Vec8 &h)
{
    double Q, Qnew, hh, bvirt;
    double enow = e[0];
    double fnow = f[0];
    int eindex = 0;
    int findex = 0;

    if ((fnow > enow) == (fnow > -enow))
    {
        Q = enow;
        enow = e[++eindex];
    }
    else
    {
        Q = fnow;
        fnow = f[++findex];
    }

    int hindex = 0;
    if (eindex < elen && findex < flen)
    {
        if ((fnow > enow) == (fnow > -enow))
        {
            Qnew = enow + Q;
            hh = Q - (Qnew - enow);
            enow = e[++eindex];
        }
        else
        {
            Qnew = fnow + Q;
            hh = Q - (Qnew - fnow);
            fnow = f[++findex];
        }
        Q = Qnew;
        if (hh != 0)
        {
            h[hindex++] = hh;
        }

        while (eindex < elen && findex < flen)
        {
            if ((fnow > enow) == (fnow > -enow))
            {
                Qnew = Q + enow;
                bvirt = Qnew - Q;
                hh = Q - (Qnew - bvirt) + (enow - bvirt);
                enow = e[++eindex];
            }
            else
            {
                Qnew = Q + fnow;
                bvirt = Qnew - Q;
                hh = Q - (Qnew - bvirt) + (fnow - bvirt);
                fnow = f[++findex];
            }
            Q = Qnew;
            if (hh != 0)
            {
                h[hindex++] = hh;
            }
        }
    }

    while (eindex < elen)
    {
        Qnew = Q + enow;
        bvirt = Qnew - Q;
        hh = Q - (Qnew - bvirt) + (enow - bvirt);
        enow = e[++eindex];
        Q = Qnew;
        if (hh != 0)
        {
            h[hindex++] = hh;
        }
    }

    while (findex < flen)
    {
        Qnew = Q + fnow;
        bvirt = Qnew - Q;
        hh = Q - (Qnew - bvirt) + (fnow - bvirt);
        fnow = f[++findex];
        Q = Qnew;
        if (hh != 0)
        {
            h[hindex++] = hh;
        }
    }

    if (Q != 0 || hindex == 0)
    {
        h[hindex++] = Q;
    }
    return hindex;
}


// 自适应高精度叉积计算
double orient2dadapt(double ax, double ay, double bx, double by, double cx, double cy, double detsum)
{
    double acxtail, acytail, bcxtail, bcytail;
    double bvirt, c, ahi, alo, bhi, blo, _i, _j, _0, s1, s0, t1, t0, u3;

    const double acx = ax - cx;
    const double bcx = bx - cx;
    const double acy = ay - cy;
    const double bcy = by - cy;

    s1 = acx * bcy;
    c = splitter * acx;
    ahi = c - (c - acx);
    alo = acx - ahi;
    c = splitter * bcy;
    bhi = c - (c - bcy);
    blo = bcy - bhi;
    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
    t1 = acy * bcx;
    c = splitter * acy;
    ahi = c - (c - acy);
    alo = acy - ahi;
    c = splitter * bcx;
    bhi = c - (c - bcx);
    blo = bcx - bhi;
    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
    _i = s0 - t0;
    bvirt = s0 - _i;
    Vec4 B = vec4();
    B[0] = s0 - (_i + bvirt) + (bvirt - t0);
    _j = s1 + _i;
    bvirt = _j - s1;
    _0 = s1 - (_j - bvirt) + (_i - bvirt);
    _i = _0 - t1;
    bvirt = _0 - _i;
    B[1] = _0 - (_i + bvirt) + (bvirt - t1);
    u3 = _j + _i;
    bvirt = u3 - _j;
    B[2] = _j - (u3 - bvirt) + (_i - bvirt);
    B[3] = u3;

    double det = estimate(4, B);
    double errbound = ccwerrboundB * detsum;
    if (det >= errbound || -det >= errbound)
    {
        return det;
    }

    bvirt = ax - acx;
    acxtail = ax - (acx + bvirt) + (bvirt - cx);
    bvirt = bx - bcx;
    bcxtail = bx - (bcx + bvirt) + (bvirt - cx);
    bvirt = ay - acy;
    acytail = ay - (acy + bvirt) + (bvirt - cy);
    bvirt = by - bcy;
    bcytail = by - (bcy + bvirt) + (bvirt - cy);

    if (acxtail == 0 && acytail == 0 && bcxtail == 0 && bcytail == 0)
    {
        return det;
    }

    errbound = ccwerrboundC * detsum + resulterrbound * fabs(det);
    det += acx * bcytail + bcy * acxtail - (acy * bcxtail + bcx * acxtail);
    if (det >= errbound || -det >= errbound)
        return det;

    s1 = acxtail * bcy;
    c = splitter * acxtail;
    ahi = c - (c - acxtail);
    alo = acxtail - ahi;
    c = splitter * bcy;
    bhi = c - (c - bcy);
    blo = bcy - bhi;
    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
    t1 = acytail * bcx;
    c = splitter * acytail;
    ahi = c - (c - acytail);
    alo = acytail - ahi;
    c = splitter * bcx;
    bhi = c - (c - bcx);
    blo = bcx - bhi;
    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
    _i = s0 - t0;
    bvirt = s0 - _i;
    Vec4 u = vec4();
    u[0] = s0 - (_i + bvirt) + (bvirt - t0);
    _j = s1 + _i;
    bvirt = _j - s1;
    _0 = s1 - (_j - bvirt) + (_i - bvirt);
    _i = _0 - t1;
    bvirt = _0 - _i;
    u[1] = _0 - (_i + bvirt) + (bvirt - t1);
    u3 = _j + _i;
    bvirt = u3 - _j;
    u[2] = _j - (u3 - bvirt) + (_i - bvirt);
    u[3] = u3;
    Vec8 C1 = vec8();
    int C1len = sum(4, B, 4, u, C1);

    s1 = acx * bcytail;
    c = splitter * acx;
    ahi = c - (c - acx);
    alo = acx - ahi;
    c = splitter * bcytail;
    bhi = c - (c - bcytail);
    blo = bcytail - bhi;
    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
    t1 = acy * bcxtail;
    c = splitter * acy;
    ahi = c - (c - acy);
    alo = acy - ahi;
    c = splitter * bcxtail;
    bhi = c - (c - bcxtail);
    blo = bcxtail - bhi;
    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
    _i = s0 - t0;
    bvirt = s0 - _i;
    u[0] = s0 - (_i + bvirt) + (bvirt - t0);
    _j = s1 + _i;
    bvirt = _j - s1;
    _0 = s1 - (_j - bvirt) + (_i - bvirt);
    _i = _0 - t1;
    bvirt = _0 - _i;
    u[1] = _0 - (_i + bvirt) + (bvirt - t1);
    u3 = _j + _i;
    bvirt = u3 - _j;
    u[2] = _j - (u3 - bvirt) + (_i - bvirt);
    u[3] = u3;
    Vec12 C2 = vec12();
    int C2len = sum(C1len, C1, 4, u, C2);

    s1 = acxtail * bcytail;
    c = splitter * acxtail;
    ahi = c - (c - acxtail);
    alo = acxtail - ahi;
    c = splitter * bcytail;
    bhi = c - (c - bcytail);
    blo = bcytail - bhi;
    s0 = alo * blo - (s1 - ahi * bhi - alo * bhi - ahi * blo);
    t1 = acytail * bcxtail;
    c = splitter * acytail;
    ahi = c - (c - acytail);
    alo = acytail - ahi;
    c = splitter * bcxtail;
    bhi = c - (c - bcxtail);
    blo = bcxtail - bhi;
    t0 = alo * blo - (t1 - ahi * bhi - alo * bhi - ahi * blo);
    _i = s0 - t0;
    bvirt = s0 - _i;
    u[0] = s0 - (_i + bvirt) + (bvirt - t0);
    _j = s1 + _i;
    bvirt = _j - s1;
    _0 = s1 - (_j - bvirt) + (_i - bvirt);
    _i = _0 - t1;
    bvirt = _0 - _i;
    u[1] = _0 - (_i + bvirt) + (bvirt - t1);
    u3 = _j + _i;
    bvirt = u3 - _j;
    u[2] = _j - (u3 - bvirt) + (_i - bvirt);
    u[3] = u3;
    Vec16 D = vec16();
    int Dlen = sum(C2len, C2, 4, u, D);

    return D[Dlen - 1];
}

double orient2d(double ax, double ay, double bx, double by, double cx, double cy)
{
    const double detleft = (ay - cy) * (bx - cx);
    const double detright = (ax - cx) * (by - cy);
    const double det = detleft - detright;

    if (detleft == 0 || detright == 0 || (detleft > 0) != (detright > 0))
    {
        return det;
    }

    const double detsum = fabs(detleft + detright);
    if (fabs(det) >= ccwerrboundA * detsum)
    {
        return det;
    }

    return -orient2dadapt(ax, ay, bx, by, cx, cy, detsum);
}

// 判断点是否在多边形内
bool pointInPolygon(const IsoBandPoint &p, const std::vector<IsoBandPolygon> &targetPolygon)
{
    int i, ii;
    int k = 0;
    double f, u1, v1, u2, v2;
    IsoBandPoint currentP, nextP;

    double x = p.x;
    double y = p.y;

    int numContours = targetPolygon.size();
    for (i = 0; i < numContours; i++)
    {
        if (targetPolygon[i].bbox.has_value())
        {
            const auto &bbox = targetPolygon[i].bbox.value();
            const auto &west = bbox[0];
            const auto &south = bbox[1];
            const auto &east = bbox[2];
            const auto &north = bbox[3];
            if (p.x <= west || p.x >= east ||
                p.y <= south || p.y >= north)
            {
                return false;
            }
        }

        ii = 0;
        const auto &contour = targetPolygon[i].coords;
        int contourLen = contour.size() - 1;

        currentP = contour[0];
        if (currentP.x != contour[contourLen].x || currentP.y != contour[contourLen].y)
        {
            throw std::runtime_error("First and last coordinates in a ring must be the same");
        }

        u1 = currentP.x - x;
        v1 = currentP.y - y;

        for (; ii < contourLen; ii++)
        {
            nextP = contour[ii + 1];

            u2 = nextP.x - x;
            v2 = nextP.y - y;
            if (v1 == 0 && v2 == 0)
            {
                if ((u2 <= 0 && u1 >= 0) || (u1 <= 0 && u2 >= 0))
                {
                    return false; // 点在边界上，返回false
                }
            }
            else if ((v2 >= 0 && v1 <= 0) || (v2 <= 0 && v1 >= 0))
            {
                f = orient2d(u1, u2, v1, v2, 0, 0);
                if (f == 0)
                {
                    return false; // 点在边界上，返回false
                }
                if ((f > 0 && v2 > 0 && v1 <= 0) || (f < 0 && v2 <= 0 && v1 > 0))
                {
                    k++;
                }
            }
            currentP = nextP;
            v1 = v2;
            u1 = u2;
        }
    }

    return (k % 2 != 0);
}

bool isInside(std::vector<IsoBandPoint> &points, const IsoBandPolygon &targetPolygon)
{
    for (size_t i = 0; i < points.size(); i++)
    {
        if (!pointInPolygon(points[i], {targetPolygon}))
        {
            return false;
        }
    }
    return true;
}

IsoBandPolygon smoothPolygon(const IsoBandPolygon &poly, int iterations)
{
    IsoBandPolygon tempOutput;
    for (int y = 0; y < iterations; y++)
    {
        std::vector<IsoBandPoint> outCoords;
        const auto &coords = y == 0 ? poly.coords : tempOutput.coords;
        for (auto k = 1; k < coords.size(); k++)
        {
            auto p0x = coords[k - 1].x;
            auto p0y = coords[k - 1].y;
            auto p1x = coords[k].x;
            auto p1y = coords[k].y;
            outCoords.push_back(IsoBandPoint(0.75 * p0x + 0.25 * p1x,
                                             0.75 * p0y + 0.25 * p1y));
            outCoords.push_back(IsoBandPoint(0.25 * p0x + 0.75 * p1x,
                                             0.25 * p0y + 0.75 * p1y));
        }
        outCoords.push_back(outCoords[0]);

        tempOutput.coords = std::move(outCoords);
    }

    return tempOutput;
}

emscripten::val gen_point(const double &x, const double &y)
{
    std::vector<double> data;
    data.resize(2);
    data[0] = x;
    data[1] = y;
    auto view = emscripten::typed_memory_view(data.size(), data.data());
    emscripten::val jsTypedArray = emscripten::val::global("Float64Array").new_(view);
    return jsTypedArray;
}

emscripten::val gen_typeArray(const std::vector<double> pointCoords)
{
    auto view = emscripten::typed_memory_view(pointCoords.size(), pointCoords.data());
    return emscripten::val::global("Float64Array").new_(view);
}

emscripten::val gen_coords(const std::vector<IsoBandPoint> &bandCoords, bool flipY,
                           double west, double south, double north, double scaleX, double scaleY)
{
    auto pointNum = bandCoords.size();

    std::vector<double> pointCoords;
    pointCoords.resize(pointNum * 2);
    for (auto i = 0; i < pointNum; i++)
    {
        double x = bandCoords[i].x;
        double y = bandCoords[i].y;
        x = west + x * scaleX;
        if (flipY)
            y = north - y * scaleY;
        else
            y = south + y * scaleY;

        pointCoords[2 * i] = x;
        pointCoords[2 * i + 1] = y;
    }

    return gen_typeArray(pointCoords);
}

emscripten::val gen_polygon(const std::vector<IsoBandPolygon> &bandRing, bool flipY,
                            double west, double south, double north, double scaleX, double scaleY, int iterations)
{
    auto ring = emscripten::val::array();
    for (auto k = 0; k < bandRing.size(); k++)
    {
        if (iterations > 0)
        {
            auto poly = smoothPolygon(bandRing[k], iterations);
            auto coordinates = gen_coords(poly.coords, flipY,
                                          west, south, north, scaleX, scaleY);
            ring.set(k, coordinates);
        }
        else
        {
            auto coordinates = gen_coords(bandRing[k].coords, flipY,
                                          west, south, north, scaleX, scaleY);
            ring.set(k, coordinates);
        }
    }
    return ring;
};